﻿#' @TODO 尝试不同的聚类方法与距离，进行聚类分析，并输出KM曲线结果
#' @title # 一致性聚类
#' @description 遍历大部分距离与方法，调用子函数`KM_curve_CC`绘制KM曲线，自动判断最佳聚类，输出结果文件到`output_dir`
#' @param exp  表达谱文件
#' @param clinical  生存信息文件需要有sample、time、status做列名,对应样本名，生存时间，生存状态。
#' @param seed 设置种子
#' @param k.max 尝试最大聚类数
#' @param output_dir 结果输出路径，需要以`"/"`为结尾
#' 
#' @return 
#' #### 聚类结果*list*
#'  - 第一层：不同聚类方法下的结果
#'  - 第二层： 
#'         - fit_list，所有的fit结果，结果是 *list*;
#'         - data_use_list，多个 列名为`sample` `cluster` `status` `time` 的数据框 ;
#'         - 最佳聚类k值，*list*
#' @return 在output 目录下生成 结果文件夹
#' 
#' @usage
#'     cluster_res <- CC_cluster(
#'        exp = exp[intrest_gene, ],
#'        clinical = clinical,
#'        seed = seed,
#'        k.max = k.max,
#'        output_dir = output_dir
#'     )
#' @export 
#' 
#' @author *WYK*
#' 
CC_cluster <- function(exp = NULL, clinical = NULL, seed = 1110, k.max = 5, output_dir = "./") {
  tmp_wd <- getwd()

  wd <- output_dir
  var_name <- deparse(substitute(exp))

  library(ConsensusClusterPlus)
  library(tidyverse)
  library(survival)
  library(survminer)
  library(gridExtra)
  # library(crayon)  
  clinical <- clinical %>% filter(time > 0,!is.na(status))
  exp <- na.omit(exp)
  common_sample <- intersect(colnames(exp),clinical$sample)
  exp <- exp[,common_sample]
  rownames(clinical) <- NULL
  clinical <- clinical %>% column_to_rownames('sample') %>% .[common_sample,] %>% rownames_to_column('sample')

  if (!dir.exists(paste0(wd, "/output/CC_cluster"))) {
    dir.create(paste0(wd, "/output/CC_cluster"),recursive = T)
  }

  title <- sprintf("CC_%s", var_name)

  getOptK <- function(res, minCls = 2, maxCls = k.max) {
    # 最佳分类数
    Kvec <- minCls:maxCls
    x1 <- 0.1
    x2 <- 0.9 # threshold defining the intermediate sub-interval
    PAC <- rep(NA, length(Kvec))
    names(PAC) <- paste("K=", Kvec, sep = "") # from 2 to maxK
    for (i in Kvec) {
      M <- res[[i]][["consensusMatrix"]]
      Fn <- ecdf(M[lower.tri(M)])
      PAC[i - 1] <- Fn(x2) - Fn(x1)
    }
    optK <- Kvec[which.min(PAC)]
    return(optK)
  }

  km_res_list <- list()
  j <- 1
  # , "minkowski","binary","maximum","canberra"
  for (distance in c("euclidean", "spearman", "pearson")) {
    if (distance == "euclidean") {
      clusteralgs <- c("pam", "hc", "km")
    } else {
      clusteralgs <- c("pam", "hc")
    }

    for (clusteralg in clusteralgs) {
      ccwd <- paste0(wd, "/output/CC_cluster/CC_", clusteralg, "_", distance)

      # print(paste0("\nRuning now: ", clusteralg, " + ", distance))
      message(stringr::str_glue('\n\nRuning now: {crayon::blue(clusteralg)} + {crayon::blue(distance)}'))

      if (!dir.exists(ccwd)) {
        dir.create(ccwd,recursive = T)
      }

      setwd(ccwd)

      res <- ConsensusClusterPlus(
        d = as.matrix(exp), # 提供的需要聚类的数据矩阵，其中列是样本，行是features，可以是基因表达矩阵。
        maxK = k.max, # 聚类结果中分类的最大数目，必须是整数。
        reps = 1000, # 重抽样的次数
        pItem = 0.8, # 样品的抽样比例，如 pItem=0.8 表示采用重抽样方案对样本的80%抽样，经过多次采样，找到稳定可靠的亚组分类。
        pFeature = 1, # Feature的抽样比例
        clusterAlg = clusteralg, # 使用的聚类算法，“hc”用于层次聚类，“pam”用于PAM(Partioning Around Medoids)算法，“km”用于K-Means算法，也可以自定义函数。
        innerLinkage = "ward.D2", 
        # the agglomeration method to be used. This should be (an unambiguous abbreviation of) one of
        # "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA),
        # "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC).
        # "single" 最短距离法 "complete" 最长距离法 "median" 中间距离法(=WPGMC）
        # "average" 类平均法（=UPGMA) "centroid" 重心法（=UPGMC) "ward" 离差平方和 "Mcquitty" 相似分析法 (=WPGMA)
        finalLinkage = "ward.D2",
        distance = distance, # 计算距离的方法，有pearson、spearman、euclidean、binary、maximum、canberra、minkowski。# nolint
        ml = NULL,
        tmyPal = NULL, # 可以指定一致性矩阵使用的颜色，默认使用白-蓝色
        seed = seed,
        plot = "pdf",
        writeTable = T, # 若为TRUE，则将一致性矩阵、ICL、log输出到CSV文件
        weightsItem = NULL, # 样品抽样时的权重
        weightsFeature = NULL, # Feature抽样时的权重
        verbose = F,
        corUse = "everything" # everything：遇到缺失数据时，相关系数的计算结果将被设为missing;
        # complete.obs：行删除;pairwise.complete.obs：成对删除，pairwisedeletion # nolint
      )

      best.k <- getOptK(res)

      km_res <- KM_curve_CC(
        res = res,
        clinical = clinical,
        var_name = var_name,
        clusteralg,
        distance = distance,
        best.k = best.k,
        k.max = k.max
      )
      
      if(j == 1){
        km_res_list[[1]] <- km_res
        names(km_res_list)[j] <- stringr::str_c(distance, "_", clusteralg)
      }else {
        km_res_list <- rlist::list.append(km_res_list, km_res)
        names(km_res_list)[j] <- stringr::str_c(distance, "_", clusteralg)
      }
      j <- j + 1

    }
  }
  setwd(tmp_wd)

  return(km_res_list)
}
#' `CC_cluster`子函数，用于绘制KM曲线与计算p值等
KM_curve_CC <- function(res = res, clinical = clinical, var_name = var_name, clusteralg = clusteralg,
                        distance = distance, best.k = best.k, k.max = k.max) {
  i <- 2
  p_list <- list()
  fit_list <- list()
  data_use_list <- list()
  best.k_list <- list()
  repeat {
    data_use <- res[[i]]$consensusClass %>%
      as.data.frame(., make.names = T) %>%
      rownames_to_column(., var = "sample") %>%
      dplyr::rename("cluster" = ".") %>%
      inner_join(., clinical)

    fit <- surv_fit(formula = Surv(time, status) ~ cluster, data = data_use)
    p <- ggsurvplot(fit,
      data = data_use,
      surv.median.line = "hv",
      palette = "lancet",
      ggtheme = theme_bw(12),
      pval = T,
      pval.size = 6,
      risk.table = T,
      pval.method = T,
      tables.theme = theme_bw(12),
      xlab = "Time",
      risk.table.fontsize = 3.6,
      legend = "top"
      # conf.int = F,
      # risk.table.y.text.col = F, # 使用颜色代替Y轴文字
      # risk.table.y.text = F # Y轴不使用文字注释
    )

    p2 <- ggpubr::ggarrange(
      p$plot+labs(x = '')+theme(plot.margin = unit(c(0.2,.2,-.15,0.2),'cm')),
      p$table+theme(plot.margin = unit(c(-.15,0.2,0.2,0.2),'cm')),
      ncol = 1,
      align = "v",
      heights = c(0.88 - i * 0.051, 0.14 + i * 0.025)
    )

    p_list[[i - 1]] <- p2

    if (i == best.k) {
      p_list[[k.max]] <- p2
    }

    if ((i - 1) == 1) {
      fit_list[[1]] <- fit
      names(fit_list)[i - 1] <- str_c("k", i)

      data_use_list[[1]] <- data_use
      names(data_use_list)[i - 1] <- str_c(clusteralg, "_", distance, "_", "k", i)

      best.k_list[[1]] <- str_c("k", best.k)
      names(best.k_list)[i - 1] <- str_c("k", best.k)
    } else {
      fit_list <- rlist::list.append(fit_list, fit)
      names(fit_list)[i - 1] <- str_c("k", i)

      data_use_list <- rlist::list.append(data_use_list, data_use)
      names(data_use_list)[i - 1] <- str_c(clusteralg, "_", distance, "_", "k", i)

      best.k_list <- rlist::list.append(best.k_list, str_c("k", best.k))
      names(best.k_list)[i - 1] <- str_c("k", best.k)
    }

    i <- i + 1

    if (i > k.max) {
      km_res <- list('fit_list' = fit_list, 'data_use_list' = data_use_list, 'best.k_list' = best.k_list)
      break
    }
  }

  cowplot::ggsave2(
    filename = sprintf("CC_%s_%s_%s_KM.pdf", clusteralg, distance, '_'),
    plot = gridExtra::marrangeGrob(p_list, nrow = 1, ncol = 1, top = ""),
    width = 6, height = 6
  )

  return(km_res)
}
#' 汇总`CC_cluster`结果，包含最佳p值、k、最优k值与KM曲线对应的p值
#' @export  
#' 2get_cc_res_table 输出结果
#' > cc_res_table
#' # A tibble: 28 × 6
#'    method        clusteralg distance  k     besk_k sur_p.val
#'    <chr>         <chr>      <chr>     <chr> <chr>      <dbl>
#'  1 euclidean_pam pam        euclidean k2    k4         0.208
#'  2 euclidean_pam pam        euclidean k3    k4         0.491
#'  3 euclidean_pam pam        euclidean k4    k4         0.613
get_cc_res_table <- function(cluster_res = cluster_res) {
    cc_res_table <- map_dfr(
        names(cluster_res),
        function(x) {
            tibble(
                # 方法与距离汇总
                method = x,
                # 方法汇总
                clusteralg = str_split(cluster_res[[x]][["data_use_list"]] %>% names(),
                    pattern = "_", simplify = T
                )[, 1],
                # 距离汇总
                distance = str_split(cluster_res[[x]][["data_use_list"]] %>% names(),
                    pattern = "_", simplify = T
                )[, 2],
                # k值汇总
                k = str_split(cluster_res[[x]][["data_use_list"]] %>% names(),
                    pattern = "_", simplify = T
                )[, 3],
                # 对应的最优k值汇总
                besk_k = cluster_res[[x]][["best.k_list"]] %>% names(),
                # 各个k值下对应的生存曲线p值汇总
                sur_p.val = map_dbl(cluster_res[[x]][["data_use_list"]],
                 function(y) {
                    surv_diff <- survdiff(Surv(time, status) ~ cluster, data = y)
                    p.val <- 1 - pchisq(surv_diff$chisq, length(surv_diff$n) - 1)
                })
            )
        }
    )
    return(cc_res_table)
}



#' @title  筛选得到最优k值
#' @description 
#' 在满足聚类后KM曲线p值条件下，筛选得到最优k值;  
#' 如果没有最优k值满足p值，则查找满足p的非最优k值;  
#' 如果都没有，说明聚类不成功，主干测试失败，屏幕打印失败信息，程序继续循环.
#' @param cc_res_table  函数`get_cc_res_table`结果
#' @param output_dir 差异基因使用聚类方法，结果记录文件路径
#' @param CC_res_sur_p 筛选聚类结果时，KM使用p值
#' @return 单行 *tibble*,筛选得到的k值与p值
#' @export  
#' @author *WYK*
#'
get_cc_res_filted <- function(cc_res_table = cc_res_table,CC_res_sur_p = CC_res_sur_p, output_dir = output_dir) {
    n_best_k <- cc_res_table %>%
        filter(k == besk_k) %>%
        filter(sur_p.val < CC_res_sur_p) %>%
        nrow()

    if (n_best_k != 0) {
        cc_res_filted <- cc_res_table %>%
            filter(k == besk_k) %>%
            filter(sur_p.val < CC_res_sur_p) %>%
            arrange(sur_p.val) %>%
            .[1, ]

        message(sprintf(
            "cluster's result in best.k is %s, method now: %s, surplot_p.val is %.4f",
            cc_res_filted$besk_k,
            cc_res_filted$method,
            cc_res_filted$sur_p.val
        )) 
    } else {
        message(sprintf("Cluster's results in best.k get no p.val < %s, Use k's result\n",CC_res_sur_p))
        n_k <- cc_res_table %>%
            filter(sur_p.val < CC_res_sur_p) %>%
            nrow()

        if (n_k != 0) {
          cc_res_filted <- cc_res_table %>%
            filter(sur_p.val < CC_res_sur_p) %>%
            arrange(sur_p.val) %>%
            .[1, ]

          message(sprintf(
            "cluster's result in k is %s, method now: %s, surplot_p.val is %.5f",
            cc_res_filted$k,
            cc_res_filted$method,
            cc_res_filted$sur_p.val
          ))
        } else {
          print(sprintf("All Cluster's results get no p.val < %s",CC_res_sur_p))
          cc_res_filted <- tibble()
        }
    }
    return(cc_res_filted)
}

# > cc_res_filted
# # A tibble: 1 × 6
#   method      clusteralg distance k     besk_k sur_p.val
#   <chr>       <chr>      <chr>    <chr> <chr>      <dbl>
# 1 pearson_pam pam        pearson  k2    k2         0.197


# source("/Pub/Users/wangyk/Project_wangyk/Codelib_YK/some_cancers/major_test2.r")


#' @TODO 聚类差异分析
#' @title  聚类差异分析
#' @description  调用多个子函数  
#' 构建基因list表达谱，先进行聚类分析，后续进行差异分析，然后单因素cox回归，lasso回归分析，构建模型.
#' 
#' @param exp 表达谱数据,是一个变量名
#' @param clinical 临床信息表数据，要求 sample status time做列名。
#' @details 临床数据的样本名，需要和表达谱中样本名一致，或者绝大部分都相同
#' @param intrest_gene 聚类使用基因集
#' @param k.max 一致性聚类最大聚类数 
#' @param CC_res_sur_p 一致性聚类后，KM曲线显著性阈值 
#' @param direction 差异基因取并集时，表征方向 
#' @param log2FC_value 差异基因阈值，绝对值 
#' @param DEG_adjusted_pvalue 差异基因阈值 校正后p值 
#' @param DEG_pvalue 差异基因阈值，p值。
#'   > 注：如果不是`NULL`，优先使用p值,如果是`NULL`，则使用校正后p值。默认`NULL`，以`NULL`为主
#' @param pval_univ_cox 单因素回归P值 
#' @param output_dir 结果输出路径，会在该路径下生成一个output文件夹，存放所有结果，需要以`"/"`结尾
#' @param seed 设置种子,默认1110
#' @param saveplot 是否保存图片，lasso结果图 与 训练集KM、ROC
#' @param surtime_unit 临床生存信息所用单位12或365
#' @export 
#' @return  *list*
#'  - 元素1：*list*，`cluster_res`，聚类分析全部结果
#'  - 元素2：*tibble*，`cc_res_table`，各聚类方法下的KM p值与最优k值
#'  - 元素3：*list*，`deg_all`，全部的差异分析结果，聚类结果是几类，就做了几次差异分析  
#'    > 举例：
#'    > - k=2, 1 vs 2; 2 vs 1
#'    > - k=3, 1 vs 2&3; 2 vs 1&3; 3 vs 1&2
#'    > - ...
#'  - 元素4：deg_all_res
#'  - 元素4：*data.frame*，`univcox_res`，单因素cox回归分析结果
#'  - 元素5：*list*，`lasso_res`，子函数lasso分析结果
#'  - 元素6：*list*，`training`，训练集结果，包含km roc曲线以及结果data.frame
#' @details 
#'   - 如果在`output`文件夹下只有 聚类结果 那么说明，聚类结果不支持差异基因分析
#'   - 如果在`output`文件夹下只有 聚类结果 + CC_res_used_for_DEGs.txt，说明，聚类结果支持差异分析，但是差异分析后，不支持单因素cox回归分析 
#'   - 如果在`output`文件夹下只有  聚类结果 + CC_res_used_for_DEGs.txt + lasso_res 说明，聚类后差异分析，单因素cox回归分析后，可以lasso，但是lasso后，不支持建模
#'   - 如果文`output`文件夹下有四个文件，聚类结果、 CC_res_used_for_DEGs.txt、lasso、训练集KM&ROC曲线，那么在条件C_res_used_for_DEGs下，聚类成功，支持后续分析，建模成功   
#' @usage 
#' major_test <- CC_DEGs(exp = brca_exp_02, clinical = BRCA_sur,
#'                    intrest_gene = inf_gene,
#'                    k.max = 5,
#'                    CC_res_sur_p = 0.2, direction = "up",
#'                    log2FC_value = 0.1, DEG_adjusted_pvalue = 0.1, DEG_pvalue = 0.1,
#'                    pval_univ_cox = 0.1, surtime_unit = 365,
#'                    output_dir = "/pub/users/innertech/wyk/Project_output/GAP454/",
#'                    saveplot = T,
#'                    seed = 1110)
#' 
#' @author *WYK*
#' 
CC_DEGs <- function(exp = NULL, clinical = NULL,
                    intrest_gene = NULL, k.max = 5,
                    CC_res_sur_p = 0.05, direction = "up",
                    log2FC_value = 0.585, DEG_adjusted_pvalue = 0.05, DEG_pvalue = NULL,
                    pval_univ_cox = 0.05, surtime_unit = c(12, 365),
                    output_dir = "./",
                    saveplot = F,
                    seed = 1110) {
  if (!dir.exists(output_dir)) {
    dir.create(output_dir, recursive = T)
  }

  var_name <- deparse(substitute(exp))

  cluster_res <- CC_cluster(
    exp = exp[intrest_gene, ],
    clinical = clinical,
    seed = seed,
    k.max = k.max,
    output_dir = output_dir
  )

  cc_res_table <- get_cc_res_table(cluster_res = cluster_res)

  cc_res_filted <- get_cc_res_filted(cc_res_table = cc_res_table, CC_res_sur_p = CC_res_sur_p)

  if (nrow(cc_res_filted) != 0) {
    write_tsv(
      x = cc_res_filted,
      file = paste0(output_dir, "/output/CC_res_used_for_DEGs.txt"),
      quote_escape = "none"
    )

    if (cc_res_filted$k == cc_res_filted$besk_k) {
      message(sprintf(
        "When method is %s, we get best.k %s,the CC_sur_p.val is %.4f.",
        cc_res_filted$method,
        cc_res_filted$besk_k,
        cc_res_filted$sur_p.val
      ))
    } else {
      message(sprintf(
        "We can't get best.k, because all sur_p.val is below %s, so when method is %s, we get k is %s,the CC_sur_p.val is %.4f.",
        CC_res_sur_p,
        cc_res_filted$method,
        cc_res_filted$k,
        cc_res_filted$sur_p.val
      ))
    }

    # 找到最佳k值下的分组的引索
    index <- list(
      cc_res_filted$method,
      str_c(cc_res_filted$clusteralg, "_", cc_res_filted$distance, "_", cc_res_filted$k)
    )

    # 找到最佳k值下的分组信息
    sur_cc_data <- cluster_res[[(index[[1]][1])]][["data_use_list"]][[(index[[2]][1])]]

    # 找到最佳k值下，聚类结果。聚成几类，循环几次，聚成两类，循环两次，取并集，不影响结果；聚成多类，循环多次，取并集。
    cc_group_index <- sur_cc_data$cluster %>%
      unique() %>%
      length()

    # 用于差异分析的分组结果，是一个*list*
    group_list <- lapply(as.list(1:cc_group_index), function(x) {
      group_for_degs <- data.frame(
        sample = sur_cc_data$sample,
        group = ifelse(sur_cc_data$cluster == unique(sur_cc_data$cluster)[x], "EG", "CG")
      )
    })

    deg_all <- lapply(as.list(1:cc_group_index), function(x) {
      tmp <- DEGs(exp = exp, group_infor = group_list[[x]], var_name = var_name) # 计算差异基因
      tmp2 <- model_data(
        degs = tmp$degs,
        log2FC_value = log2FC_value,
        DEG_adjusted_pvalue = DEG_adjusted_pvalue,
        DEG_pvalue = DEG_pvalue,
        exp_1 = tmp$exp_1
      ) # 按阈值筛选差异基因
      return(tmp2$deg2)
    })

    deg_all_res <- map_dfr(deg_all, function(x) {
      if (direction == "up") {
        x %>%
          filter(log2FC > 0)
      } else if (direction == "down") {
        x %>%
          filter(log2FC < 0)
      } else {
        print("Param direction must be \"up\" or \"down\".")
      }
    })

    deg_gene_res_chara <- deg_all_res %>%
      select(gene) %>%
      pull() %>%
      unique()

    if (length(deg_gene_res_chara) > 1) {
      major_res <- univCox2KM_ROC(
        exp = exp[deg_gene_res_chara, ],
        clinical = clinical,
        pval_univ_cox = pval_univ_cox,
        seed = seed,
        saveplot = saveplot,
        output_dir = output_dir,
        surtime_unit = surtime_unit,
        var_name = var_name
      )

      if (length(major_res) > 1) {
        tmp3 <- list(
          cluster_res = cluster_res,
          cc_res_table = cc_res_table,
          deg_all = deg_all,
          deg_all_res = deg_all_res,
          univcox_res = major_res$univcox_res,
          lasso_res = major_res$lasso_res,
          training = major_res$training
        )
      } else {
        tmp3 <- list()
      }
    }
  } else {
    tmp3 <- list()
  }

}


